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ABSTRACT 

A numerical model was developed depicting the wind-driven 
circulation in an ocean basin containing an island. This 
linear, barotropic, filtered model was utilized to test and 
evaluate the "Hole Relaxation" technique (Allen 1954) in 
preparation for future comparative studies with a free sur- 
face model and later incorporation in a multi-level world- 
ocean model with islands. Eight different data cases were 
studied to evaluate the model’s ability to properly treat a 
variety of island siz 
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I . INTRODUCTION 



For decades, men driven by the reality of the world’s 
oceans and their geometric configurations have strived to 
model the physical boundary characteristics and the processes 
which govern the ocean-atmosphere realm. These men have 
recognized the future needs of the world’s forecasting media 
by developing models which have answered some old queries 
and have further related new theories to make present-day 
phenomena more easily understood. Most studies were con- 
siderably enhanced by the advent and sophistication of 
numerical techniques . 

In 1947, Sverdrup showed for the first time the major 
role of the curl of the wind stress in determining the 
general mid-ocean distribution of the meridional current. 
Later, Stommel (1948) demonstrated the decisive role of the 
latitudinal change in the coriolis parameter in the forma- 
tion of intensive boundary currents along the eastern shores 
of continents. Stommel also used the linearized vorticity 
equation with linear bottom friction. In 1950, Munk, utiliz- 
ing lateral frictional dissipation, directed his studies at 
more detailed features of linear viscous boundary currents, 
all the while restricting these studies to simple domains of 
integration (i.e., rectangle and triangle) and by zonal wind 
fields. All of these models used a vorticity equation con- 
sidered on a Beta plane. 
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'‘"Further studies, now advanced by numerical techniques, 
were undertaken. Gates in 1968 showed the development of 
the large-scale wind-driven circulation using the primitive 
equations of motions for a homogeneous ocean model charac- 
terized by a free upper surface. This time-dependent model 
permitted not only easy observation of the individual tran- 
sients’ role in the circulation development, but also exten- 
sive examination of the role of divergence in the large-scale 
circulation by the use of the free surface upper boundary 
condition. Bryan in 1968 and later in 1969 with Manabe, by 
means of a joint atmosphere-ocean numerical model with a 
rigid upper boundary in the ocean, investigated the role of 
the ocean circulation in producing realistic features of the 
climate of the atmosphere. 

Also in 1969, Takano studied the barotropic ocean circu- 
lation in a doubly connected world ocean using the vorticity 
equation (non-s tationary) and the numerical method of "hole 
relaxation." Kamenkovich et_ _al. (1969), in an independent 
effort, further displayed successful simulation of actual 
shoreline profiles by numerical methods in a complete circu- 
lation of the world ocean (stationary). In 1971, Alexander 
extended Gates 1 work to a baroclinic free surface model depict 
ing a doubly connected domain, with islands in the western 
boundary region. 

The primary purpose of the following study is to test a 
numerical technique for calculating the large-scale wind- 
driven circulation in an ocean basin which contains an island. 
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It is not the purpose of this model to physically depict any 
particular ocean or seasonal wind pattern; rather, the goal 
is to simplify the circulation problem in order to study and 
evaluate the numerical technique at the sacrifice of precise 
world simulation. 

The numerical solution of a barotropic ocean with an 
island has been achieved in both linear and non-linear models. 
The present work is totally a linear effort utilizing the 
vorticity equation which filters the external high frequency 
gravity waves. This permits the use of a longer time step 
(=14 hours) than in the free surface (non-f il t er ed) models. 

The application of the n Hole Relaxation” technique (Allen 
1954) with respect to the island boundaries (interior 
boundaries to the ocean model) is employed with a major goal 
to test and evaluate the techniques on this simple model for 
future development. Furthermore, the results of this 
evaluation will be used in calculating the barotropic mode 
in more sophisticated world ocean models which can ill- 
afford the luxury of short time steps as required in the 
free surface models. 

The present study did not attempt to argue in favor of 
a filtered model rather than a free surface model for the 
calculation of the barotropic mean currents. Its goal is to 
simply test a technique whereby the merits of the two 
numerical approaches (free surface vs. rigid surface) may be 
compared . 
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II. THE MATHEMATICAL STATEMENT OF THE PROBLEM 



A. VORTICITY EQUATION DERIVATION 

The problem is to obtain a solution to the wind driven 
ocean circulation in a rectangular ocean basin of uniform 
depth containing an island. 

The linearized equations of horizontal motion (1) and 
(2) together with the continuity equation (3) for a homo- 
geneous fluid are as follows: 



3u , .1 3p ,n2 • 1 3T : 

vr - fv + — — IT*- = AV U + — — — T — 

9t p 3x p 3z 



3v , , , 1 3p .„2 , 1 

tt— + f u + — — = AV v + — 

3t p 3y 



3t 



p 3z 
o 



3u 3jv 3w _ n 

3x 3y 3z 



(1) 

( 2 ) 

(3) 



Friction forces are crudely represented by a lateral 

eddy diffusion of momentum (AV 2 u,AV 2 v) and a vertical eddy 

stress (x ,T ). All other symbols have their usual meaning 
x y 

Cross differentiation of (1) and (2) to eliminate pressure 
results in the following vorticity equation: 



3 ,3v 3u. c , . ,3u 3v, _ 

^ +vB +f( aX + ' 

A (-i — V 2 v - V z u)+ -i 
v 3x 3y P Q 3z 3x 

„ 3f 

where 3 - -r — = constant 

3y 



3x 

37 



, (4) 



By vertically integrating the continuity equation (3) from 
the bottom of the ocean (z = -D) to the mean sea level (z=0) , 
we have 



3u , 3v_ , w(o) - w(-D) _ 0 

3x 3y D 



(5) 
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where Cu,v) is the vertically integrated velocity defined by 



- 1 

u 58 r- / u dz , 

-D 

and 

_ 1 0 

V = zr / V dz . (6) 

-D 



The boundary conditions on w(-D) and w(o) are 

w(-D) = 0 

w ( o ) = 0 (7) 

The first boundary condition in (7) is kinematic and 
applies at the flat ocean floor , The second condition is 
the approximation which removes external gravity waves from 
the model and thereby allows longer time steps to be 
utilized. The inaccuracy of the solution due to the balance 
approximation of requiring the vertical velocity to be zero 
at the sea surface can be examined by comparison of the 
results of this present calculation with similar calculations 
using a barotropic primitive equation model with a free sur- 
face as used for example by Gates in 1967, As a result of 
(7) , (5) becomes 



2“ + iJL = o 

3x 3y 



(8) 



Consequently, a streamf unction i^(x,y) is defined for the 
vertically integrated velocity as follows: 



u ' 3y 
v = H 

9x 



(9) 



The vertical integration of the vorticity equation (4) and 
the use of (8) and (9) results in 
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1 9l v 
T P q D 1l 3x 3y 



3t 

' ) z = 0" ( 3x 



3T x 

3f>z = .D^ • 



Cio) 



The boundary conditions on the stress at the ocean floor 
(z - -D) and at the sea surface (z = 0) are as follows: 
At z = -D 



while at z = 0 




T 

X 



- F cos 

b 

0 , 



(ID 



( 12 ) 



where F is the amplitude of the zonal components of the wind 
stress and b is the north-south extent of the ocean basin. 
This model in the islandless case was a time dependent repre- 
sentation of Munk’s model (1950). Therefore, the results 
can easily be compared with Munk’s analytic solution as a 
check on the accuracy of the numerical solution. 

Thus, there is no bottom stress, and the surface stress 
corresponds to a pattern of westerly winds in the northern 
half of the basin, and to easterly winds in the southern 
half of the basin. In addition, it is important to have a 
stress that has a non-zero curl since only the curl of the 
wind stress actually does the forcing in the vorticity equa- 
tion. Using (11) and (12), the vorticity equation takes the 
final form, 

- -off - ^Db sln <?> . (13) 
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B. BASIN COASTLINE AND ISLAND BOUNDARY CONDITIONS 



The boundary conditions along the island and basin coast- 
lines (mainland) extremities used in the solution of (13) are 
those of zero normal flow and of zero slip. 

The condition of no normal flow (v^ = 0) was implemented 
by requiring that the s t r eamf unc t ion be equal to zero on the 
ocean periphery and that the same function be a constant 
function of time on the island’s perimeter. These conditions 
are written as : 

^ = 0 (mainland boundary) 

iJj = c(t) (island boundary) (14) 

The value of the constant c(t) determines the net flow 
(volume transport) between the mainland and the island and 
is calculated by the "hole relaxation” procedure (Allen 1954) 
as a part of the solution. 

The conditions of zero slip and no normal flow are 
implemented in the calculation of the lateral friction force 
term on the right hand side of (13). These conditions are 
written as : 

u, v = 0 (all island and mainland 

boundaries). (15) 

Details of the formulation of these boundary conditions are 
contained in section three where the finite difference equa- 
tions are given. 

C. PHYSICAL CHARACTERISTICS OF THE NUMERICAL MODEL 

The model presented here is a rectangular basin of 

, L = 8000 km; of breadth, b = 4000 km; and of constant 



length 
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depth. Depth - 2 km, (See Figure 1.) The total grid size is 
41 x 21 in the horizontal while the computational grid size 
is 40 x 20, The distance between adjacent grid points is 
200 km in both x and y. The value chosen for Beta makes the 
grid centered at 32.5°N. Figure (1) also displays a portion 
of the 861 nodal points and shows how an arbitrary island 
may appear. The arbitrarily sized, rectangularly shaped 
island may be placed at any selected location throughout the 
grid. The model is not strictly designed to simulate any 
particular ocean or ocean currents system, but rather is 
utilized to evaluate the solution techniques to be tested 
and secondarily to evaluate the evolution of the physical 
flow patterns around the island with respect to time. 

D. EXPERIMENTS CONDUCTED 

Initially, an islandless ocean model was studied. This 
model utilized the same equation (13) as its basis of solu- 
tion. Again, this model was analogous to the solution of 
Munk’s model (1950) when the numerical solution achieved 
steady-state conditions. All characteristics were as de- 
scribed above with the exception of an island. A duration 
of 210 days was selected as the run time for this phase in 
order to be able to compare all possible "island” integra- 
tions with this islandless case. 

Next, an island-ocean basin phase was introduced by 
adapting the boundary conditions and restrictions of the 
hole relaxation technique to be tested (Allen 1954). A 



15 



{ HH°N-s»» j6 ®a } 3 P n »!IM 



Q 

O 

in 



in 

<N 

ro 



Q 

n 



•ooo#ooo#c#ooo#o#ooooooooo#ooo#o#o#o#ooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•oooooooooooooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•oooooooooooooooooooooooogoooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooogooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•oooooooooooooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo© 

ooooooooooooooooooooooooooooooooooooooooo 

•oooooooooooooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

•oooooooooooooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooooooooooooooo 

•ooooooooooooooooooooooooooooooooooooooo# 

ooooooooooooooooooooooooooooooooooooooooo 

to#o«oto«o#o«o®o#oeo»o«o«o«o«o#o#otototo« 

OOOOOOOOOOOOOOOOOOOOeOQOOOOOOOGOQOOOOOQOO 

•ooooooooo# dogogogogo- • ' 0 o#o#oto« oooonoo# 
O 0 O#O 0 O 0 O 0 O 0 L OOOOO •GOOOGOOOGOOOOOOOO# o#o 
*ooooooooo#f># :>#c •oooogoooooooqooogooooooo 

ooooooc© • • e • © • © ooo ©o© • ©o#o 
•ooo#oo*#ci#q©o#g#gooo • ogo-oogocogoo© • • 

OOOOOOO# ©• • ^©O© • O 3 O O W < © © • • • D 

• 0#0#0#0OG0C#G©<; ©OOOOi 3#0#0#G# 

000©0©0© '• © e • o • * • O O • • s © « «) 
•OOOOOOOOO#^# • "•'©. ©G© • '• • '© ?• .OOOOOO* • 

0*0©O#0#G0 r >0C • ROCOCO ^©OOGOCOCOf •0#G#C ‘#0#0 

•OOOOO# D#G#0#0#0#0#0#G#0#0#0#000#0# 2#0#0# 

o#o0ofo0ot^0n0o«o#o#o«o#o0o0 gooooo }#*:>•■ ;oo 

•0#Q©0#0#0#0OG#O o#co?#o#.# #o#o# 

OOOOOOOOOOOOGOOO #0#0#U#0# 0# ->• OOO 

•OOO0O0O0O0^#G#O 0#0#G#0#0#0#0#0# 

OOOOOOOOOOCi© iOO# •O#O#OOO#0OO#QOO 

•0#0«0#Q#0#0#0#G OOOOOOOOOOOf DOOO 

o©o©o©o©G© r ©o© ©o©o©o©?©o© • • ) 

•o©o©o©o©o©g©o©o 0©0©0©0©0©Gr©G©£© 

o©o©o©o©o©o©c ©o© ©o©o©o©o©o©o©o©o 

•o*o©o©o©o©g©0©o GO D©C • ■ • ©n©0©0© 

0©0©G©0© 3>©G© i>©G©0©C©Q©0©0© ' • • • ^© GOOOCOO 

•0©0©0©G©0©0© ©0©n©0©0©0© j© ‘ • ?©G • ?©C©G© GO 

0©0©0©0©0©C©0©0©*. ©o©o©o©c© :oooc ©G©0©0©0©0 

•0©0©0©0©G©G©0©0©0©G© :-©Q©0©G©0© '©G©G© j©G© 

o©o©o©o© o©©©o©o© • • © • © • • ©o 
©o©o©o©o©^© •••••©•©••••••• 

0©0©0©G©0©0©'3© • 0©G©!j©0©O©O©O©C©O©3©G© ©o 

•o©o©o©o©o©o©o©o©o©o©o©o©o©o©o©o©o©o©o©o© 



m 



CD 



to 



w I — ) -f-' 

a <2 *3 



rO £ 



0) 

c *p 



T3 



§ 



-P cv 

CO OG 



<D 



•C ^ 



CO 



P, c 



ro 



CN 



I “'’Ic 01 ( M»P!M P! J 9 






16 



time plot of the total kinetic energy showed that 210 days 
was adequate to reach a statistical steady state, and to 
give an accurate comparison to the steady-state characteris- 
tics of the islandless phase results. 
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III. FINITE DIFFERENCE EQUATIONS 



Individual terms on the right hand side of (13) are ex- 
pressed in finite difference form in the following paragraphs. 
The grid points in the annulus are discussed first in each 
case, then follows a discussion of the grid points on the 
island perimeter. Island corner points and the points on the 
perimeter adjacent to the corners are not discussed. Readers 
are referred to the documented computer program in Appendix A 
for further details. It should also be clarified that all 
terms on the right hand side of (13) were calculated after 

multiplying the entire equation by the area represented by a 

2 

grid point, namely (Ax) . The reason for this action was 



that 


the time 
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of change 
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as seen 
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endix 
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It should 


be po 


int ed 


out 


tha 


t a forcing 


f unc t ion 


of 


F1 = 0 



was assumed (Allen 1954) in the interior nodes of the island. 
Finally, since it was assumed that Ax = Ay throughout the 
grid, the term "Ax" was utilized throughout the subsequent 
discussion and all sections of this thesis. 

A. THE BETA TERM 

In the annular locations while maintaining boundary con- 
dition (14) , the beta term was expressed as 

„ 3<i> * 2 „ ,^1+1 " ^I-1 N A ,, 

-g Ax ^ -g ( 2 ) Ax (16) 
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If the beta terms on the north and south perimeters were 
considered to be zero due to boundary condition (14), namely, 
zero normal flow, then the boundary conditions of constant 
values of i p on the perimeter of the island and zero values of 
ip on the mainland, inherent in (14), had to be maintained on 
the east and west island perimeter nodal points. The con- 
straint that the integrated sum of all beta terms along any 
latitudinal plane or axis was to be zero, guided the formula- 
tion of the finite difference form of the beta term on the 
east-west perimeters. This integral constraint which was a 
property of the differential expression for also had to 

o X 

hold true in the finite difference analogues, as long as the 
meridional velocity component on the island perimeter was 
evaluated by means of an uncentered difference approximation. 



Therefore, when the 


beta term 


was 


calculated 


on the 


island 


perimeter as describ 


ed above, 


the 


uncentered 


f ini te 


difference 


forms used were: 













-3 || Ax 2 =^-3 (^ - iJ’j.-l) Ax , (17) 



for the western perimeter, and 

-3 f| Ax 2 4- -3 0J> I+1 - 1 > T ) Ax , (18) 

for the eastern perimeter. 



B. CURL OF WIND STRESS TERM 

In both the case of the annular location and 
island perimeter locations, the wind stress term 
in the same manner, namely: 



that of the 
was defined 
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TTF a 2 

p Db Ax sin 
r o 






- b Ax sin (Y > » 
r o 


(19) 


where 




Y 


IT (J-l ) Ax 
b 


(20) 



C* CURL OF FRICTION FORCE TERM 



Boundary conditions utilizing zero slip and zero normal 
flow (14) were again selected to enable dissipation to occur 
at the boundaries through the linear curl of the friction 
force term in the model* 

The friction force term was expressed as follows: 



[AV^V 2 *)]^- [AV 2 C1i>j - tA(f^ 2 v - | 7 V 2 u)] I jJ . (21) 

The finite difference form in the open oceanic areas 
away from the boundaries was based on (21) as 



V 2 v T +V 2 v 



t (AV 2 C) IjJ Ax 2 ]=^AAx[ ( ^ 



I , J-l 



)-( 



V 2 v -f* V 2 v 

I-l.J V I-1,J-1 

O ) 



V 2 u t ,+V 2 u. 



» , v V 2 u + V^u 

-(— Liill), , (22) 



where 



V 2 v 



V -4- V +v -hv — 4 V 

- / I+l.J+2 I+l.J 1+2, J+l I,J+1 I+1,J+1 X 

I,J 1 Ax 2 ' ’ 



(22a) 

with a similar expression for V 2 u T , and where u T and 

J. y J 1 J J 

v . are defined in terms of \p T as, 

1 9 J 1 j J 

1 ^I — l J — l"^^I J — 1 ^1 — 1 T j 

U I,J = Ax ^ ( 1,J 2 } } ~ ( 2 ~ ~^) ] , (22b) 

and 

1 $ T t T — 1 ^ T — 1 T^"^T — 1 J— 1 

V I,J = A^~ [ ( ? 2 1 * 2 ~ " )] * (22C) 
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On the island perimeter, the same basic finite difference 
form (21) applied except that the following alterations 
existed : 

i) on the northern and southern perimeters. 



These conditions preserved boundary condition (14) with 
respect to no normal flow. 

Having made use of (23) and (24) along with the zero slip 
condition, the following finite difference forms resulted 
for use in (22) on the island perimeter: 
i) southern perimeter of the island 




(23) 



and ii) on the eastern and western perimeters, 




(24) 




(25) 



ii) northern perimeter of the island 




“ U I, J+2~ U I+1, J+2^ ’ 



(26) 



iii) western perimeter of the island 




and iv) eastern perimeter of the island 



(22)4f ; [-4v 




(28) 
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This complex perimeter calculation has been clearly 
illustrated in Figure 2. 

Suppose the point (I,J) in Figure 2 was assumed to be 
located on the western perimeter of the island. The rela- 
tionship, for example, 



cfevM T , - ( 



I,J 



V 2 v +V 2 v, V 2 v ? +V 2 v 1 

^ ") - ( ~) 



(29) 



demonstrated how the staggered grid points, numbered one 

through 12, were used to compute (*r — V 2 v) . 

0 X I j J 

The right hand side of (29) was written as the sum 

2 w t * • v T . . where w T . . T1 . was the net weighting 

A > J 

factor assigned to v at the nodal point (I + i, J + j) in 

9 2 

the computation of v at the nodal point (I,J) in the 

island perimeter. Having rewritten the right hand side of 
(29), and having used the zero slip condition, v n 
it was found that: 






-v 



I-i,J+j f 

(29)^>4v t _+4v 
i , j 



I, J+l V I , J+2~ V I-1 , J+1~ V I-1 , J~ V I , J-l 



(30) 



Readers should note that in the staggered grid depicted, 

v was located at point number one, v - was located at 
1 , J J. y J * -L 

point number two, etc., and that the terms of (30) were 
signed as a result of (29),i.e., points one and two had 
w = +4, while points nine through 12 had w = -1. 

The above served to demonstrate the derivation of (27). 
Equations (25), (26), and (28) were derived in a similar 

fashion. 
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Island Perimeter 

Key : 

• s Grid Nodal Points 

o : Staggered Nodal Points 

• = ^ 

° s u,v,Vu r vv 

Figure 2; Illustration represents a portion of 

the computational grid, which includes 
the island perimeter. The side of the 
grid which represents ocean or island 
depends on the case to be evaluated, 
i.e., on the right for the Equation 
(29) cited. The shifting and/or re- 
orientation of the island portion of 
the figure will permit north, south, 
east or west perimeter calculations. 
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D. OTHER TERMS 



In addition to the three terms on the right hand side of 
(13), the finite difference forms for the residual calcula- 
tion in the relaxation phase of the model ! s solution were to 
have been expressed. Yet prior to that discussion, a brief 
description of the time differencing procedures was formu- 
lated. 

For accuracy, a centered time difference scheme was used. 
This was the so-called "leapfrog 11 scheme which was charac- 
terized by the desirable property of being neutral and the 
undesirable property of producing a computational mode in 
time for the system solution. This resulting mode, caused 
by the presence and utilization of three time levels in the 
description of , was removed by the periodic employment 

of a Euler-Backward (Matsuno) time scheme. The Euler-Back- 
ward scheme was used initially at the first time step and 
again at every 50 time steps thereafter to prevent separation 
of the solution in time. This scheme removed the computa- 
tional mode by utilization of only two time levels for the 
calculation of -|“j~ . 

The relaxation phase of the model solution had an equa- 
tion that was to be solved which was expressed as, 

V 2 (|£) - FI = 0 (31) 

where FI was expressed by the right hand side of (13) multi- 
2 

plied by Ax . If it is assumed that the symbol \p ' was the 
first derivative of ip with respect to time, then the relaxa- 
tion residual, RESID, was expressed as 
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;• 




•* 













( 32 ) 



in both the ocean basin grid points and at the island grid 
points. 



E. RELAXATION DESCRIPTION 



1 • Open Ocean Area 

The Liebmann over-relaxation technique was utilized 

in all open ocean areas of the computational grid. This 

technique was found to have the most, successful application 

to the model in which forcing by large scale stress dominated. 

Following the standard relaxation technique, if it 

was found that in the solution of (31) the absolute value of 

RESID in (32) was greater than some small allowable error, 

1 , J 

£, a new ifi’ denoted by if;' was defined as follows: 

x , j x , J 



♦i.' E *i,j + 



(33) 



where the correction C was determined by the condition 

1 , J 

* * 

that the new residual value RESID calculated using if;’ T , 

i , J i > J 



must be a negative fraction of RESID- _ , i.e., 

X , J 

RESID* = - a RESID 

X , J X , J 



(34) 



A value of a = 0 therefore corresponds to no over- 
relaxation. Since 

a 



(35) 



is true, then it was easily found that the relationship for 

the correction C _ was 

I , J 
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I ^T.q a ] RESID 



C36) 



C 






X,j 



An optimum value for the over-relaxation coefficient, 
a, of 0.75 was selected for use in these open oceanic regions 
This value was the result of several test evaluation runs 
and was rather large due to the large scale nature of the 
forcing and resultant motion. More details as to the appli- 
cation of this technique to the open oceanic areas of the 
grid can be found in the ample commentary of the computer 
program located in Appendix A. 

2 . Island Areas 

In the island areas, the "Hole Relaxation" technique 
was used (Allen 1954). This technique assumed that the same 
Equation (31) was to be solved with the exception that the 
island points were cumulatively treated like a block where 
the residual on the island, RESIDI, was defined as 



RESIDI = 



E 



(V ip 1 - Fl I>a ) , 



(37) 



island 

T symbol represented a four grid point 

average, and the symbol T / j 1 represents the finite dif- 

island 



where the T - 



ference approximations to an integral over the island of 
the individual residual values. This total residual value, 
RESIDI, was made to approach zero in compliance with the 
hole relaxation technique in order to achieve an accurate 
solution for (31) . As a result of the averaging process 
used, RESIDI was written as the sum of the individual nodal 
point residuals, each multiplied by a weighting factor where 
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the corner perimeter nodal points were assigned a factor of 
0.25, while the lateral perimeter points were given a factor 
of 0.50, and the interior points a factor of 1.0. 

If the subscripts 1,2, 3, 4; E,W,N,S; and INT. were 
used to denote the contribution to RESIDI from the four 
corners, the east, west, north and south lateral perimeter 
points and the interior nodal points respectively, then 
since FI = 0 on the island interior, (37) was expanded as: 
RES IDI=0 . 25 [ (V 2 ^'-F1) 1 +(V 2 ^ , -F1) 2 +(V 2 ^ , F1) 3 +(V 2 ^-F1) 4 ] 

+ 0.50[Z(V 2 if>'-Fl) + Z (V 2 i|>'-F 1)+Z (V 2 ip ' -Fl)+E (V^'-Fl) ] 

,E W N S 

+1.0[ I (V 2 ^ 1 )] . (38) 

INT. 

In accordance with the hole relaxation technique 
(Allen 1954), and the boundary conditions of no normal flow 
(14), had the same value at all points on the island, i.e 

E tp f (on the island) . (39) 

l , J 

The object of the hole relaxation technique was to 

find the value of for which the absolute value of RESIDI 

was less than the error tolerance E . The procedure used 

to do this was as follows. After completing each relaxation 

sweep over all interior ocean points, RESIDI was calculated 

from (38) and its absolute value tested. If it was equal to 

or less than e , no correction was added to ip 1 , and if no 

correction had been made to any of the interior ocean points 

on that last sweep, then ij;’ and comprised the solution. 

l , J 

If the absolute value of RESIDI was larger than c , a correc 
tion. Cl, was added to and a new relaxation sweep over 
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all grid points was initiated. The correction Cl was deter 

mined in such a way that the new RESIDI would be a negative 

fraction of the current RESIDI. The correction Cl for the 

island block was determined similarly to the development of 

C_ for the ocean points as follows: 
i , J 

It was first noted that in (38), 

a. for each corner grid point, 

(V 2 ^ f ) = -2\J; T + non-island iJj* ; (40) 

n-l ,2,3,4 I,J 

b. for each group of lateral nodes not on a corner 

for example the eastern perimeter, 

E(V 2 ip’) = E (-ip ' )+E (non-island ip ’ ) , (41) 

E EE 1,J 

or 

E (V 2 ip ’) = (~\p ') (NO„) + E (non-island ip ’ ) , (42) 

E E i,J 

where NO was the number of non-corner nodal points on the 

£j 

island f s eastern perimeter. Finally 

c. for the nodal points on the interior of the 

i s land , 

O^V^j = 0 . (43) 

If the domain of the island was described by desig- 
nated I,J indices such that the minimum and maximum values 
of I and J were IMIN and IMAX, JMIN and JMAX respectively, 
then the number of nodal points on a particular island side 
were expressed in terms of these indices. The number of 
non-corner lateral nodal points in the example stated in 
(42 ) was , 
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(44) 



N0 £ ~ JMAX - (JMIN + 1) 

Similar expressions for westerly, northerly, and 
southerly lateral side descriptions were evolved for NO^, 

NO^, and NOg, respectively. 

When Equations (40), (42), (43), and (44) and their 

associated counterparts for the other perimeters were used 
to expand (38), the resulting expression for the island 
residual was , 

RESIDI = 0 . 25 [ -8^ ' + E (non-island ipl ) -E (FI) ] 

n=l , 2 , 3 , 4 I,J n n = l ,2,3,4 n 

+0 . 50 [ -t/j ' (2 IMAX-2 IMIN-2+2 JMAX-2 JMIN-2 ) 

+ E (non-island ip' ) - E (FI) ] , 

n=E,W,N,S i,J n n=E,W,N,S n 

which was simplified to 

RESIDI = -(TERM)(^ ! )+ Z (non-island jp 1 ) , (45) 

where 

TERM = IMAX - IMIN + JMAX - JMIN . (46) 

The "hole relaxation" stipulated that individually 
there was no need for restriction on the values of island 
residuals, provided only that their sum, namely RESIDI, was 
eventually brought to zero. Therefore, when the absolute 
value of (45) was greater than the tolerable error e , a new 
value of the island ip’, denoted i p' , was defined by 

\p 1 * = ij; ' + Cl , (4 7) 

where Cl was determined by requiring that the new island 

* , 
residual, RESIDI , be a negative fraction of RESIDI, that is, 

RESIDI* = -(al) RESIDI , (48) 
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where al was the appropriate over-relaxation coefficient for 
the chosen island. Since from (45) 

RE SID I = (TERM) C^'*) + Z (non-island i|> ’ ) , (49) 

by substituting (47) into (49) and using (48), the constant 
Cl was obtained to be 

CI = " ‘ t erm 0111 RESIDI (50) 

The optimum value of the over relaxation coefficient, 
OtX, was found to depend upon the size of the island selected 
for study. When the island sizes were small with respect to 
the scale of'the dominant wind stress, a value of 0.40 was 
assigned to al. If the island sizes were comparable to the 
scale of the stress, a smaller value of 0.20 was assigned. 
These time-tested selections were based upon the comparative 
need for more or less over-relaxation respectively as the 
island size was varied from small to large. 

More precisely, greater over-relaxation was neces- 
sary when the scale of the island was small in comparison 
to the scale of the wind stress and less over-relaxation was 
found necessary when the scale of the island was comparable 
to the scale of the wind stress. These results were to be 
expected since the scale of the motion, and hence the scale 
of the residuals will be determined by the scale of the wind 
stress forcing . 

Further details as to the application of this tech- 
nique in conjunction with that of the open oceanic areas 
can be found in Appendix A. 
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IV. RESULTS 



A. A RESTATEMENT OF 
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evaluate the 
The accuracy 
be evaluated 
face model . 
poration will be 
with islands . 

The reader is reminded tha 
model’s solution are already r 
multi-level baroclinic models. 

In the programming of the 
made to minimize computer time 
value of the over-relaxation c 
to be near the optimum value, 
required 28 minutes of C.P.U. 
integrate 364 time steps. Thu 
approximately 14 hours, each o 
shown in the following section 
one hour (.47 hours) of comput 
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B. GENERAL DESCRIPTION OF ISLAND DATA CASE STUDIES 

Altogether, there were eight studies considered. Of 
these studies, there were seven different island positions 
examined, six of which were in the vicinity of the western 



31 



boundary; while the other island position was in the eastern 
boundary vicinity. The remaining study showed the island- 
less case. Because of the linearized form of the equations 
and the barotropic assumptions, the primary motivation for 
these eight studies was to specifically test and evaluate the 
model techniques and arrive at some quantitative idea as to 
what to expect in a more general, non-linear, multi-level 
future application of this model without excessive computer 
cost. 

Each island case study contained a variation of three or 
more island sizes centered about a specifically chosen grid 
location. For example, in the second data study, a specific 
grid position of '’WESTERN GRID HALF-CENTERED" was chosen. 

The three distinct island sizes comprising this particular 
set were 400 kilometers square, 800 kilometers square, and 
1200 kilometers square. All three islands in this set were 
centered in the western portion of the grid so that symmet- 
rical s treamf unc t ion patterns would appear to the north and 
south of the island. Each of the other island case studies 
was similarly structured and can be referred to in Table I. 

C. EQUILIBRIUM CONSIDERATIONS AND FINDINGS 

A most important aspect of the numerical solution was to 
determine if the model had achieved a steady-state condition 
in each case study. The existence of an equilibrium condi- 
tion was accurately evaluated by means of a total kinetic 
energy (T.K.E.) calculation (Appendix A). At the end of each 
time step, this mathematical indicator projected an integral 
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picture of the steady-state characteristics of the evaluated 
island case. Specific illustrative comparisons between the 
steady-state characteristics of the islandless case and 
those of the island cases are shown in Figures 3 through 9. 
From an examination of these figures, it was clear that a 
statistical equilibrium was achieved in all cases by 210 days 
of integration. The streamline maps for the 210th day for 
each case is shown in Figures 10 through 34, while some of 
the data on each case is summarized in Table I. 

As shown in Table I, the following trends were noted: 
a) When the island was located in a region of strong currents 
(Cases 2 , 3 , 4 , 6 , 7 , 8) , the larger the size of the island, the 
lower the final value of the total kinetic energy. (Case 8-d 
was an unexplained exception to this rule.) When the island 
was located in a region of weak currents (Case 5), the stream- 
lines were easily bent around the island (Figure 23) pro- 
ducing a stronger flow in the central open ocean. Because 
the kinetic energy depends on the square of the current 
speed, this streamline bending by the island increased the 
total kinetic energy. b) Except for Case 8-d, the closer 
the island was positioned to the mainland in the westerly 
located studies, the smaller the final value of the total 
kinetic energy. c) In Case 6 and in Case 7, the placement 
of a larger island in the path of a strong current reduced 
the final total kinetic energy value. d) The comparative 
examination of Case 3 and Case 4 pointed out that the same 
final T.K.E. values were achieved by the symmetric placement 
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of equivalently sized islands south and north of the central 
latitude. In Case 4^b, an asymmetrically placed island 
showed a different T.K.E. value which still conformed, how- 
ever, to those conditions stipulated above in (a) and (c) . 
e) Finally, as might be anticipated, all island cases of 
comparatively small size had final T.K.E. values similar to 
that achieved in Case 1, the islandless case. 

As an additional guide in the above interpretation, the 
average kinetic energy was calculated from the T.K.E. by 
dividing by the number of non-island values which entered 
into the total. The results were unchanged; namely, that the 
larger the island, the smaller the average kinetic energy in 
all cases except the Eastern Island Case 5 and Case 8-b, 
where a large average kinetic energy appeared with the larger 
island . 

The contents of Table I display among other things the 

value of the streamf unction on the island (il>. - and maxi- 

r island 

mum streamf unction value in the ocean basin (to ). These 

max . 

values were significant in that they presented two quantita- 
tive measures of the current’s volume transport per unit 
depth. Because of the zero normal flow boundary condition 
(14) , displayed the actual value of the volume trans- 

port per unit depth between the island and the mainland. 

The reader is especially cautioned to note (in Case 8) that 
because of recirculation between the westerly boundary and 
the island, this value of ^ s i an( j should be interpreted only 
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as a measure of the net transport between an island and any 
mainland . 

The if) value accurately pointed out the maximum volume 

transport per unit depth to be found in the ocean basin. This 

value always was found in the area immediately adjacent to 

the western boundary of the mainland. 

As shown in Table I, in each data case study: a) The 

values of ip. were smaller in magnitude than those of 

island & 

ip . b) When the size of the island was relatively larger, 

the magnitudes of ip . - , and Jp were both smaller (Cases 

island max. v 

2, 3 and 4). c) When the island was positioned closer to 

the southern mainland boundary in the westerly located 

studies, the values of it. - , and those of \p were both 

island max . 

smaller. In the easterly located studies (Case 5), the 

ip . - , values were nearly all the same, but the ip values, 

island ' r max. 

found near the western boundary of the ocean basin, consis- 
tently increased with the size of the island. This \p 
increase seemed to be related to the increased gradient of 
the s t reamf unc t ion near the central latitudinal axis of the 
basin; however, a further, more satisfactory explanation has 
not yet been found. d) In both Cases 3 and 4, the corres- 
ponding values of il. , , and \p were identical, again 

demonstrating the symmetrical characteristics of the linear 
solution. e) Finally, in Cases 6 and 7, smaller transports 
were recorded around the islands which were fully extended 
in the paths of strong currents, while the \p values 

found near the western ocean basin boundary were seen to be 
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significantly smaller only in the case where the island was 
fully extended north to south. 

Many features of the different solutions can only be 
seen in the streamf unction isoline maps which are displayed 
in Figures 10 through 34, 

Figure 10 showed the islandless case in which the numeri- 
cal solution compared favorably with the analytic solution of 
Munk . Figures 11 to 13 showed that when an island was placed 
in the center of an ocean gyre, the solution was not changed 
appreciably until the island size was about the same size or 
larger than the gyre itself. 

One of the interesting features found in Cases 3, 4 and 

5 (Figures 14 through 23) was that in this linear model, when 

an island was placed in the southern (northern) part of the 

basin, the latitude of maximum meridional transport shifted 

northward (southward). This shift was imperceptible in Case 

8 (Figures 30 through 34) in which the island had only a 

small north-south extent. There was some form of western 

intensification evident on the eastern shore of the island 

in every case study. Yet the maximum s treamf unction values 

found on the eastern side of the island were always smaller 

in magnitude than the ^ ' values depicted in Table I, 

max , 

which were always found adjacent to the western boundary of 
the basin. The island's eastern shore intensification was 
most evident in Case 6 (Figures 24 through 26) in which 
there was an anticyclonic gyre in the center of the basin, 
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near 



the northern 



and two western intensified cyclonic gyres 
and southern ends of the fully extended island. 
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TABLE I: Data Case Studies used and results obtained over a 210 day integration 
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Duration oF Time [day*] 

Figure 3’ Total Kinetic Energy for Case 2 for island sizes 
a) 400 km on a side; b) 800 kn on a side; and 
c) 1200 km on a side, compared to Case 1, the 
islandless case 
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Duration oT Time [days] 

Figure 5: Sane as Figure 3 except for Case 4, island size 

a) 400 km on a side; b) 800 x 1200 km^; 
c) 800 km on a side; d) 1200 km on a side 
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Figure 6: Same as Figure 3 except for Case 
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Figure 8: Same as Figure 3 except for Case 7, island sizes 

a) 800 x 1600 kin 2; b) 800 x 1400 km 2 ; and 
c) 800 x 1200 kn 2 
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210 DAYS 
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Figure 1£: Same as Figure 10 except for Case 2b 
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Figure 13: Same as Figure 10 except for Case 2c 
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50 



Figure 14: Same as Figure 10 except for Case 3a 
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Figure 15s Same as Figure 10 except for Case 3b 



210 DAYS 
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Figure 16: Same as Figure 10 except 
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10 except for Case 4b 
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Figure 21 s Same as Figure 10 except for Case 5a 
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Figure 22: Sane as Figure 10 except for Case 5b 
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»ame as figure 
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Figure 27: Same as Figure 10 except for Case 7a 
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Figure 29: Sane as Figure 10 except for Case 7c 
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Figure 30: Same ss Figure 10 except for Case 8a 
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Figure 32: Sane as Figure 10 except for Case 8c 
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Figure 34: Sane as Figure 10 except for Case 8e 



v. 



CONCLUSIONS 



This numerical model, depicting the wind-driven circula- 
tion in an ocean basin containing an island, was successful 
in portraying the boundary conditions and relaxation tech- 
niques to be evaluated. Furthermore, a firm numerical basis 
has been established for a future non-linear modification and 
still later comparative study with a free surface model. The 
study of the combined relaxation technique employed in this 
model has established a much needed data resource which dis- 
plays the usefulness of the "Hole Relaxation 11 technique 
(Allen 1954) for the problem of the flow around an island. 
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, , _ _ APPENDIX A 

// EXEC FORTCLGP, REGION. G0=160K 
/ /FORT .SYS IN DD DS NAME = S SP3 ( CONTUR ) , D I SP = SHR 
// DD * 

THE MAIN PROGRAM INITIATES THE SOLUTION OF THE 
VORTICITY EQUATION BY ZEROING ALL DIMENSIONED 
VARIABLES. THE SUBRCUTI NE S, NAMELY CALF AND RELAX, THEN 
TAKE CONTROL AMD PERFORM THEIR RESPECTIVE FUNCTIONS. 
AFTERWARDS , THE MAIN PROGRAM TAKES CONTROL AND PER- 
FORMS THE TIME-STEPPING PHASE OF THE SOLUTION . 



D VAR 



MEANS DIMENSIONED VARIABLES 



************* 



***** DPSIDI IS THE VALUE OF DPSIDT ON THE PERIMETER 
****** OF THE ISLAND . THIS VALUE IS CONSTANT . ***** 

COMMON/ I SLAND/ I MIN, IMAX , JMI N, JMAX , DPS ID I 

COMMON /D VAR /PS I M 1, P SI , F 1 , U , V , DE LSQU , DE LS QV ,DPSIDT, 

IRES ID, PSTEMP 

DIMENSION PSI Ml ( 41 ,21 ) , PS 1(41, 21), FI (41, 21), U( 42, 22), 
IDELSQUI 40, 20 ) ,DELSQV( 40, 20) , DPSIDT (41 ,21 ) ,RESID( 41 ,21) 
2 , PSTEMP (4 1,21) ,V (42,22 ), DAT AH 41, 21 ) 

READ( 5,5009) I MIN, I MAX , JMI N , JMAX 
5009 FORMAT ( 414 ) 

RE AL*8 TITLE1 ( 12 ) 

REAL *4 CL ( 25) 

LOG I CAL* 1 LTG( 3) /.TRUE., .TRUE. , .FALSE./ 

DATA CL/-2. 4, -2. 2, -2. 3,-1. 8, -1.6, -1.4, -1.2, -1.0, -0.8, 
1~ 0.6, — 0.4,— 0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 
21.8,2.0,2.2,2.4/ 



START THE ZEROING PROCESS. 

DPSIDI =0.0 
DO 3000 J= 1, 21 
DO 3000 1=1,41 
PSIMlt I , J) =0.0 
PSI ( I , J ) =0 . 0 
FI (I , J) =0.0 
R E S I D ( I , J)=0.0 
PST EMP ( I , J ) =0.0 
3000 DPSIDT ( I ,J) =0.0 
DO 4000 J= 1 , 20 
DO 4000 1=1,40 
DELSQV( I ,J) =0.0 
40 00 DEL SQU ( I , J ) =0. 0 
DO 5 COO J = 1 , 22 
DO 5000 1=1,42 
U( I , J ) =0 .0 
5000 V(I,J)=0.0 

LET’S GO TO THE MATSUNO SCHEME TO GET THIS PROJECT 
GOING . 

NPR I NT = 50 

DE LT AX = 2 . 0 * 10. 0**7 
DELTAT= 5.0 *10. 0**4 
TIMEST=210 .0*86400. 0 
TI ME=0. 0 

WRITE(6, 103) TIME 

LEAPFROG TIME SCHEME ******************************** 



THE LEAPFROG SCHEME IS NEUTRAL IN CHARACTER BUT THE 
PRESENCE AND UTILIZATION OF THREE TIME LEVELS IN THE 
DESCRIPTION OF THE FIRST DERIVATIVE OF PSI WITH 
RESPECT TO TIME , PRODUCES A COMPUTATIONAL MODE IN 
TIME. THIS MODE HAS TO BE AND IS REMOVED BY THE E'JLER- 
BACKKAP D SCHEME. THIS SCHEME USES TWO LEVELS IN TIME 
IN ITS DESCRIPTION OF THE FIRST TIME DERIVATIVE OF 
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PS I .THEREFORE, THERE IS NO COMPUTATIONAL MODE IN TIME 
AND A MORE ACCURATE SOLUTION FOR THE PSI FIELD AT 
TIME LEVEL • N+ 1 ' IS ATTAINED . 

DO 3100 K = 1 , 1 0 00 
IF(K.EQ.l) GO TO 2000 
L=M0D( K , NPRINT ) 

IF (L.EQ.0) GO TO 2000 
CALL CALF 
CALL RELAX 
DO 6000 J=2 , 20 
DO 6000 1=2, AO 
TEMP=PS I ( I , J ) 

PS I ( I , J ) = PSIM1 ( I , J )+2 .0 * DELTAT * DPSIDT(I,J) 

6000 PSIM1( I , J) = TEMP 

PRINT THE PSI FIELD IN TAEULAR FORM ****************** 
GO TO 2800 



AT THIS TIME , PSI IS KNOWN AT TWO CONSEQUTIVE TIME 
LEVELS AND IS STORED IN PSI AND PSIM1 FIELDS 
RESPECTIVELY. BY MEANS OF THE ' TEMP ' ST ATEM.ENT , THE 
TWO TIME LEVELS ARE ADVANCED BY ONE INCREMENT OF TIME 
, I . E . , DELTAT = 5. 0*10. 0**4 (SEC.) RESPECTIVELY . 



THE EULER-BACKWARD OR MATSUNO TIME SCHEME ; NOTE: PSI 
AT TIME LEVEL 'N' IS SAVED IN PSTEMP FIFLD FOR 
LATER USE IN SUBSEQUENT BACKWARD IMPLICIT STEPS.. 



2000 WRITE(6,103) TIME 
2100 DO 2200 J=2 , 20 

DO 2200 1=2,40 

PS IM1 ( I , J)=PSI ( I, J ) 
2200 PSTEMP ( I , J) =PSI ( I ,J) 



BY MEANS OF THE ABOVE LOOP ( #2200 ), PS I AT TIME LEVEL N 
IS PUT INTO THREE FI ELDS , NAMELY PSI, PSTEMP AND PSIM1. 
ALL FIELDS ARE AT THE SAME TIME LEVEL IN ORDER TO 
PRESERVE LINEAR COMPUTATIONAL STABILITY IN THE 
FRICTION TERM OF THE FORCING FUNCTION, FI. 

COMMENCE THE FORWARD TIME STEP PHASE OF THE MATSUNO 
SCHEME . 

NOTE: ALL TERMS OF CALF AND RELAX ARE OUTPUTED AT TIME 
LEVEL ' N ' . 

CALL CALF 
SCALE=10 .0** ( - 2 ) 

DO 2910 J=2 , 20 
DO 2910 1=2,40 
2910 RESID( I, J)=F1( I ,J)*SCALE 

WRI TE ( 6 , 108) ( (RESID(I,J),I=1,33),J=1,21) 

108 FORMAT! • O' , 'RESID FIELD = FI FIELD * SCALE ',/ ,21 ( 33F4 

1 . 2 ,// ) ) 

CALL RELAX 

DPSIDT IS NOW AT TIME LEVEL 'N' . ******************* 

DO 2500 J =2 , 2 0 
DO 2500 1=2,40 

PSI(I,J)= PS I ( I , J ) + DELTAT * DPSIDT(I,J) 

2500 PSIMlt I , J ) =P SI ( I , J) 

PRESENTLY, THE FIELDS PSI,PSIM1 AND PSTEMP HAVE 
CONTAINED IN THEM, PSI VALUES AT TIME LEVELS • N+l' 
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INTERMEDIATE, 'N+l' INTERMEDIATE AND *N' RESPECTIVELY. 



NOW AS A PART OF THE MATSUNO TIME SCHEME ; A SIMULATED 
(IMPLICIT) BACKWARD TIME STEP IS EMPLOYED . THE 
BACKWARD STEP IS USED AS A CALIBRATION STEP TO REFINE 
THE PSI FIELD AT TIME LEVEL 'N+l' . NOW COMMENCE THE 

BACKWARD TIME STEP PHASE OF THE MATSUNO SCHEME . 

NOTE: ALL TERMS OF CALF CALCULATED AS A FUNCTION OF 
PSI AT INTERMEDIATE TIME LEVEL 'N+l' . 



CALL 

CALL 



CALF 

RELAX 



DPSIDT AT TIME 
NOW OUTPUT PSI 



LEVEL 'N+l' (INTERMEDIATE) OUTPUTED. * 
AT TIME LEVEL 'N+l'. ****************** 



I , J ) 
J) + 



DO 2700 J =2 » 2 0 
DO 2700 1=2,40 

PS I Ml ( I, J) = PSTEMP( 

2700 PSI ( I , J ) =PSTEMP ( I 
2800 TIME= DELTAT*< 

WRITE(6,103) TI ME, K 
SCALE= 10.0** (-8) 

DO 2900 J=2, 20 
DO 2900 1=2,40 
R ES I D ( I , J ) = PSI ( I ,J) * 
WRITE(6,104) ((RESIDU 



DELTAT 



DPS IDT ( I, J) 



2900 

103 

104 



SCALE 

, J ) , I = 1 , 33 ) 



•RESID FIELD = PSI 



, J = l,21 ) 

16 ,// ) 

FIELD * SCALE' ,/,21 (33^4 



FORMAT!' ' ,T63»'TIME=' ,F10.1 
FORMAT’! ' ( 

1 . 2 ,//)) 

LL=MOD ( K ,3 ) 

TDAY=TIME/86400.0 
WRITE (6 , 103 ) TDAY , LL 

LOOP IS FOR THE CONVERSION OF DATA TO A FORM MORE 
USEABLE IN THE COMPUTER SUBROUTINE , ' CGNTUR ' 

J= 1 * 21 



2950 



9970 

9999 

99 80 
3100 



DO 2950 
JS=22-J 
DO 2950 1=1,41 
DATA1( I »J)=RESID( I , J) 

IF ( (TDAY. GT. 200.0) .AND. (LL .EQ.O ) ) GO TO 9970 
GO TO 9980 

RE AD ( 5 , 9999 ) ( T I TL E 1 ( J ) , J = 1 , 12 ) 

FORMAT ( 6 A8 ) 

CALL CONTUR! DATA 1,41, 2 1,41 ,C L , 25 , TI TLE 1 ,4,08,LTG) 

IF(TIME.GE.T IMEST ) STOP 

CONTINUE 

STOP 

DEBUG 

AT 3100 

TRACE ON 

DISPLAY K 

END 



SUBROUTINE CALF 



THIS SUBROUTINE CALCULATES THE RIGHT HAND SIDE OF THE 
VORTICITY EQUATION WHICH IS MULTIPLIED BY THE SQUARE 
OF DELT AX* £ £ £ £ * ?*: ^ 'k & # ± 3$c # ^ ^ 3$c:$s x. £ # 

DVAR MEANS DIMENSIONED VARIABLES ************* 

** THE ENTIRE GRID IS DIVIDED INTO AREAS IN ORDER TO 
** PERMIT BETTER EASE FGR- CALCULATIONS IN BOTH ISLAND* 
** RELATED RESIGNS AND OCEANIC REGIONS . ************ 

COMMON/ I S L AND/ I M IN , I M AX , J M I N , J M AX , DP S I D I 
COMMON /D VAR /PS I Ml , PSI , F 1 , U , V , D EL SQU , DELSQV , DPS IDT, RE SI 
1 D, PST EMP 



oooooooooo ooooooo non o on on non on onooooonn 



DI ME NS I ON PSIM1(41,21),PSI(41,21),F1(41,21),U(42,22), 
1DELSQU ( 40,20),DELSQV(40,20),DPSIDT(41,21),RESID(41,21) 
2,PSTEMP(41,21) , V ( 42 , 22 ) 



+ PSIM1(I+1,J) 
PSIM1 < I, J) 



DEFINE INITIAL U(40,20) AND V<40,20) AS FUNCTIONS OF 
PSIM1(41,21) .REFERENCE SHOULD BE MADE TO SKETCHES OF 
GRIDS IN THE JOURNAL NOTES : ***********************- 

U(I + 1,J + 1) AND Vd+ltJ + 1) TRANSFORM U,V(40,20) 

INTERNAL GRID SYSTEM INTO EXTERNAL U,V (42,22) GRID 
SYSTEM • *****************************************:; 

DELTAX=2 .0^1 0. 0**7 
A= 10. 0**8 
DO 50 1=1,40 
DO 50 J = 1 ,20 
AA= PSIMKI + l, J + l) 

BB= PS IM1( I, J+l ) + . . . 

CC= PSIMKI + l, J+l) +PS I Ml ( I » J+ 1 ) 

DD= PSIMKI + l, J) + PS I M 1( I , J ) 

EE = 2 . 0*DELT AX 

V(I + 1,J+1) = {(AA/EE) — C BB/ EE) ) 

50 U(I + 1,J + 1) = UDD/EE) -(CC/EE)) 

NOTE: FOR ALL FRICTION FORCE TERMS ; 

CALCULATION OF FRICTION TERM USING PSIM1 

(FORWARD TIME SCHEME). *************************** 

FRICTION TERMS APE CALCULATED USING A FORWARD TIME 
SCHEME. THIS IS DONE TO PRESERVE LINEAR COMPUTATIONAL 
STABILITY. K THE LEAPFROG SCHEME WAS UTILIZED, THE 
FRICTION TERM WOULD BE UNSTABLE. ******************** 



DEFINE EXTERNAL BOUNDARY VALUES OF U AND V *********** 

DO 60 1=2,41 

U( I , 1) = -U( I, 2) 

V( I ,1) = -V (1,2) 

U( I , 22) = -U( I ,21) 

60 V( 1,22 ) = -V (1,21) 

NOTE THAT THE FOUR EXTERNAL CORNER VALUES OF U$V ARE 
NOT DEFINED BJT A VALUE OF ZERO IS ASSUMFD FOR EACH. 

IF THEY ARE TO BE DEFINED THEN U ( 4 2 , 1 ) = J ( 4 1 , 2 ) *,U ( 1 , 1 ) = 
U( 2, 2) ; U( 1, 22) =U( 2 , 21) 5 AND U ( 42 , 2 2 ) =U ( 41 , 2 1 ) ******** 



DO 70 J=2 , 21 
U(l, J)=-U(2,J) 

V( 1, J) =-V ( 2 , J ) 

U( 42 , J ) =-U ( 41 , J ) 

70 V ( 42 , J ) =-V ( 4 1 , J ) 

************** AREAS ONE AND FOUR CALCULATIONS: ****** 
************** AFEA ONE IS THE AREA WHICH DEPICTS THE * 
************** 00 FAN I C REGION SOUTHWEST , SOUTH AND **** 
************** SOUTHEAST OF THE ISLAND. AREA FOUR IS ** 
************** TH E AREA WHICH DEPICTS THE OCEANIC **** 
************** REGI ON NORTHWEST, NORTH AND NORTHEAST *** 
************** OF THE ISLAND . *********************** 

BETA TERM CALCULATION ******************************** 

BETA=1.94*10.0* 5 M-13) 

M= JMI N-l 

IF(M.LT.2) GO TO 27 
DO 25 J = 2» M 
DO 25 1=2,43 

25 Fl(I,J)=-(3ETA)*(PSI(I+l,J)-PSI(I-l,J))MDELTAX)/2.0 
27 MM= JMAX +1 

IF (MM. GT. 23) GO TO 33 
DO 30 J =MM ,20 
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DO 30 1=2,40 

30 Fl( I , J) =-( BETA)* ( PSI ( I + 1 , JJ-PS I (I-1»J) ) * (DELTAXJ/2.0 
STRESS CURL TERM CALCULATION ************************* 



THE AMOUNT OF POSITIVE VORTICITY INTO THE OCEAN FROM 
THE ATMOSPHERE IS NEGATIVE ( ANT I -CYCLON I C ). THIS 
MAKES THE OCEAN SPIN A/C LIKE A H I GH , PRDDUC I NG A 
POSITIVE STREAMFUNCT ION, PS I . ********************** 

STRESS CURL TERM CALCULATION ************************* 

33 F= I . 0 

PI=3. 1415926 
DEPTH=2. 0*10. 3**5 
B = 4. 0*10. 0**8 
M= JM I N- 1 

IFIM.LT.2) GO TO 47 
DO 45 J = 2, M 
DO 45 1=2,40 

Y=P I * ( J-l)*( 1. 0/B)*DELTAX 
STRESS= ( ( F*P I ) / ( DEPTH*B ) )*SIN( Y) 

45 FI ( I , J ) = F1 ( I, J )-(STRESS*DELTAX**2) 

47 MM = JMA X + 1 

IFCMM.GT .20) GO TO 44 
DO 40 J=MM, 20 

DO 40 1=2,40 , „ 

Y=PI*(J — 1)*{ 1 . 0/ B ) *DEL TAX 
STRESS = ( (F*PI )/ ( DEPTH* B ) )*S IN(Y } 

40 F1(I ,J)=F1( I ,J ) — ( STRESS*DELTAX**2) 

DEFINE EXTERNAL BOUNDARY VALUES OF U AND V . ********* 
44 N = I M I N + 1 

DO 62 I =N , IMAX 

U( I , J M I N + l ) =- U ( I , JMIN) 

V( I , JM I N+ 1 ) =-V ( I , JMI N) 

U ( I , J M A X ) =- U ( I , JMA X+ 1 ) 

62 V ( I , J MAX ) =- V { I , JMAX + 1 ) 

BASED ON THE (40,20) DIMENSIONED DEL SQU AND DELSQV SYS 
-TEM, THE NEXT LOOP CALCULATION IS PERFORMED . ****** 

DO 800 1=1,40 

DEFINE D EL SQU AND DELSQV. NOTE: BOTH TERMS ARE DEFINED 
WITH THE FACTOR DE L TAX-SQUAR ED CONSIDERED IN EACH 
TERM . X; £ 5k # ## # 5k st ^ ^ ^ ^ ± ^ 

DELS3U( I , J ) IS THE LA PLACI AN OF U ;DELSQV(I,J) IS THE 
LAPLAC I AN OF V. BOTH ARE DEFINED IN THE INTERIOR RE- 
GION OF THE GRID . ******************************** 

M= J M I N- 1 
DO 600 J =1 , M 

DELS3U( I ,J)=(U(I+i,J+2)+U(I+l,J)+U(I+2,J+l )+U(I,J+l)-( 
14. 0*U ( 1+1, J+l) ) ) 

600 DELSQV (I , J ) = ( V ( I + 1 , J+2 ) +V ( I + 1 , J ) +V ( I +2 , J + 1 ) +V ( I , J + 1 ) - ( 
14. 0*V( 1 + 1, J+l) ) ) 

DO 800 J = J MAX, 20 

DELSQU ( I , J) =(U ( I +1 , J+2 )+U( 1+1 , J ) +U ( 1+2, J+l ) +U ( I , J+l)-( 
1 4. 0*U ( 1 + 1, J + l) ) ) 

800 DELSQV ( I , J ) = ( V ( I +1 , J +2 ) +V ( I + 1 , J ) +V ( I +2, J + 1 ) + V ( I , J + 1 ) - ( 
14. D*V( 1 + 1 , J+l) ) ) 

CALCULATE THE CURL OF THE FRICTION FORCE , DEFINED ON 
THE GRID INTERIOR . ********************************** 

CURL OF THE FRICTION FORCE TERM. ******************* 

M= JMI N- 1 

I F ( M . LT .2) GO TO 97 
DO 95 J = 2 , M 
DO 95 1=2,40 



76 



non non ooooooooo oonooo 



AAA= ( DE LSQV ( I , J ) + DELSQV ( I , J- 1 ) ) /2.0 
BBB= ( DELSQV! 1-1 ,J) + DELSQV ( I -1 , J-l ) ) /2.0 
CCC= ( DELSQU! I , J ) + DE LSQU ( I -1 , J ) ) /2.0 
DDD= { DELSQU ( I , J-l ) + DEL SQU ( I- 1 , J- 1 ) ) /2.0 

DELVX IS THE PARTIAL OF THE LAPLAC I AN OF V WITH 
RESPECT TO X. *********************************** 

DELUY IS THE PARTIAL OF THE LAPLACI AN OF U WITH 
RESPECT TO Y « ************************************* 

DELVX= AAA - BBB 
DELUY= CCC - DDD 

95 Fl(I,J)=Fl(I,J)+(( DELVX-DELUY ) / ( DELTAX) ) *A 

97 MM= J MA X + 1 

I F ( MM . GT .20 ) GO TO 98 
DO 90 J =NM , 20 
DO 90 1=2,40 

AAA= ( DELSQV ( I , J ) + DELSQV ( I , J- 1 ) ) /2.0 
BBB= ( DELSQV! I -1 ,J) + DELSQV ( I -1 , J-l ) ) /2.0 
CCC= ( DELSQU ( I , J ) + DELSQU! I-1,J ) ) /2.0 
DDD= ( DELSQU ( I , J-l ) + DELSQU ( 1-1 , J-l ) ) /2.0 
DELVX= AAA - 3 BB 
DELUY= CCC - DDD 

90 F 1 ( I , J ) =F1 ( I , J ) + ( ( DELVX-DELUY )/ t DELTAX) )*A 

************** AREAS TWO AND THPEE CALCULATIONS: ***** 
************** AREA TWO IS THE AREA WHICH DEPICTS THE* 
************** OCEANIC REGION WEST OF THE ISLAND. **** 
************** AREA THREE IS THE AREA WHCIH DEPICTS ** 
************** THE OCEANIC REGION EAST OF THE ISLAND.* 

BETA TERM CALCULATION ******************************* 

98 BETA=1. 94*10. 0**(-13) 

DELT AX = 2 .0*10. 0**7 
NN= I M I N-l 

IF(NN.LT.2) GD TO 28 
DO 26 1=2, NN 
DO 26 J = JMI N , J MAX 

26 Fit I , J )=-(3ETA)*(PSI(I + l,J) — PS KI-l, J))*( DELTAX) /2.0 
28 MM = I MAX + 1 

I F ( MM . G T . 40 ) 30 TO 34 
DO 31 I = MM , 40 
DO 31 J = JMI N, J MAX 

31 F 1 ( I ,J)=-{BETA)* ! P SI (I +1 , J ) - PS I (1-1 ,J) ) * (DELTAX) /2.0 

STRESS CURL TERM CALCULATION : ********************** 

34 NN= I M I N - 1 

I F ( NN . L T. 2 ) GO TO 48 

DO 46 1=2, NN 

DO 46 J = JM I N, J M AX 

Y=PI*( J-l)*( 1. 0/B) *DELTAX 

STRESS= ( ( F*P I )/(DEPTH*B) )*SIN!Y) 

46 F 1 ( I , J ) =F1 ( I , J ) -( STRE$S*DELTAX**2) 

48 MM= I MA X + 1 

IF( MM . GT .40 ) GO TO 49 

DO 41 I =MM , 40 

DO 41 J=JMIN,JMAX 

Y=PI*(J-1)*(1.0/B)*DELTAX 

STRESS= (( F*P I ) / ( DE PTH* B) ) * SIN(Y) 

41 Fl(l,J)=Fl(I,J)-( STRESS*DELTAX**2) 

DEFINE EXTERNAL EOUNDARY VALUES OF U AND V . ********* 

49 MN= JM I N + 1 

DO 71 J = MN , JMA X 

U( IMIN + 1 ,J ) =-'J( I MIN, J ) 

V! IMIN + 1 , J ) =- V ( IM IN, J ) 

U ( I M A X , J ) =-U! I MAX + 1 , J ) 

71 V! IMAX, J )=-V( IMAX+l.J ) 

MX=JMAX-1 
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DO 801 J = J MI N, MX 
NN= I M I N- 1 
DO 601 I =1 , NN 

DEL SOU ( I , J ) =( U( 1 + 1 , J + 2) +U( I +1 , J)+U( 1+2 , J+l)+U ( I , J + l )- 
1 (4 . 0*U ( 1+1, J+l ) ) ) 

6 01 DEL SQ V ( I , J ) = ( V ( I +1 , J+2 ) + V ( I + 1 , J ) + V ( I + 2 , J+ 1 ) +V < I , J + 1 ) - < 
14. 0*V ( 1 + 1, J + l) ) ) 

DO 801 I = I MAX » 40 

DELSQIH I , J) = ( U( 1 + 1 , J + 2)+U( 1+1 , J) +U( 1+2 , J + l ) +U ( I , J + l >- 
1 (4.0*U( 1+1, J+l ) ) ) 

8 01 DE LSQV { I , J) = (V ( 1+1 , J+2 ) +V ( 1+1, J )+V ( I +2, J + l )+V( I , J + 1 ) - ( 
14. 0*V{ 1 + 1, J+l) ) ) 

CURL OF THE FRICTION FORCE TERM . ******************** 
NN= I MI N- 1 

IF (NN. LT.2) GO TO 92 
DO 96 1=2, NN 
DO 96 J=JMIN,JMAX 

AA A= ( DE LSOV ( I , J ) + DELS QV { I , J-l ) ) /2.0 
BBB= (DELSQV(I-1,J) + DEL SQV ( I - 1 , J- 1 ) ) /2.0 
CCC= ( DELSQU ( I , J ) + DE LSQJ ( I -1 , J ) ) /2.0 
DDD= { DELSQU( I , J-l ) + DE LS QU ( I -1 , J-l ) ) /2.0 
DELVX= AAA - BBB 
DELUY= CCC - DDD 

96 F 1 { I , J ) =F 1 { I , J) + ( ( DELVX-DELUY)/ (DELTAX) )*A 

92 NM= I MAX + 1 

IF (NM.GT.40) GO TO 93 
DO 91 I =NM , 40 
DO 91 J = JM IN, J MAX 

AA A= ( DE LSQ V( I , J ) + DE LS QV ( I , J-l ) ) /2.0 
BBB= ( DELSQVt I - 1 , J ) + DEL SQV ( I - 1 , J-l ) ) /2.0 
CCC= ( DE LSQU ( I , J ) + DELSQU ( I - 1 , J ) ) /2.0 
DDD= (DELSQUII , J-l) + DELSQU ( I -1 , J-l ) ) /2.0 
DELVX= AAA - BBB 
DE LUY = CCC - DDD 

91 FI ( I , J) =F1( I , J )+ ( ( DELVX-DELUY) / (DELTAX) )*A 

* ****** ******* AREA FIVE CALCULATIONS: *************** 
**************AREA FIVE IS THE AREA WHICH DEPICTS *** 
*************3 NLY THE PERIMETER OF THE ISLAND . ** 

** BETA TERMS FOR AREA FIVE'S NORTH AND SOUTH PERI ME-* 
** TERS ARE ZERO . ********************************* 
** WHEN J =JMI N CR JMAX AND CALCULATIONS ARE PER- ** 
** FORMED AT THE ISLAND CORNERS FOR THE BETA TERM, ** 
** THE ACTUAL FINITE DIFFERENCE c ORM IS PSIU+l) ** 
** —PS 1(1—1) DEVIDED BY 2.0, BUT BECAUSE P S I ( I + 1 ) = *** 
** PSI ( I ) =CONST ANT FOR THE SW/NW CORNERS AND PS I (I) * 
** =PSI ( I-1)=C0MSTANT FOR THE SE/NE CORNERS, THE **** 
PRESENT FORM OF FINITE DIFFERENCE CAN BE USED ** 

** VALIDLY FOR SIMPLICITY . 

** COMMENCE BETA TERM CALCULATIONS FOR AREA FIVE'S ** 
** EASTERN AND WESTERN PERIMETERS. ******************* 

BETA TERM CALCULATION ******************************** 

93 BETA=1. 94*10. 0**(-13) 

DELTAX = 2. 0*10. 0**7 
FACT = 1.0 

DO 32 J = JM I N» J MAX 

I F ( ( J. EQ.JMIN) .OR. (J.EQ.JMAX) ) FACT = 0.5 
I = I M I N 

F 1 ( I , J ) = -(BETA)*(PSI(I,J)-P SI(I-1, J) ) * ( DE LTAX) *FACT 
I =1 MAX 

32 Fill , J)=-( BETA)* (PSI ( I +1 , J ) —PS I (I, J) )*( DELTAX )*F ACT 

** COMMENCE CALCULATION OF WIND STRESS CURL TERM ON ** 
** THE NORTHERN AND SOUTHERN PERIMETERS . ************ 

N= I MI N+l 
F=l. 0 

PI=3. 1415926 
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DEPTH=2. 0*10. 3**5 

B= 4. 0* 10. 0**8 

1 1 = I MAX- I M I N 

I F { 1 1 . L E . 1 ) GO TO 11 

NX=I MAX- 1 

DO 42 I =N , NX 

J=JMI N 

Y=PI*( J-l)*( 1. 0/B) *DELTAX 
STRESS = { ( F*PI ) / ( D EPTH* B ) ) *S IN { Y) 

Fit I , J) =-( STRESS*DELTAX**2) 

J = J MAX 

Y=PI * ( J-l)* (1. 0/B)* DEL TAX 
STRESS=( ( F*P I ) / ( DE PTH*B ) )*SIN(Y) 

42 F1(I, J )=-( STRESS*DELTAX**2) 

** COMMENCE CALCULATION OF WIND STRESS CURL TERM ON ** 
** THE EASTERN AND THE WESTERN PERIMETERS . ********** 

11 DO 43 J = JMI N, J MAX 

Y = PI*(J-1)*( 1 . 0/ B ) *DEL TAX 
STRESS = ( { F* P I ) / ( DEPTH* B J )*SIN(Y ) 

I = I M I N 

FI { I i J ) = FI ( I » J )-( STRESS*DELTAX**2) 

I = 1 MAX 

43 Fl( I, J ) = F 1 { I, J )-( STRESS*DELTAX**2) 

****AREA FIVE FRICTION FORCE TERM CALCULATION ******** 
****N0TE: ALL CALCULATIONS DO NOT INCLUDE CORNER GRID* 
****POINTS . THESE CALCULATIONS WILL BE COMPLETED LAST 
****F0LL0WING CALCULATIONS ARE FOR NORTHERN AND SOUTH- 
PERIMETERS ONLY . #*£************:***** 

****TEST FOR ZERO NODAL POINT ISLAND . *************** 

I I = I MAX - I M IN 

I F ( I I . L E . 1 ) GO TO 233 

****LOOP FOR N/S PERIMETER LESS CORNERS : ************ 

N = I M I N + 1 

NX= I MAX- 1 

DO 343 I =N , NX 

I F ( I .EQ. I MIN + 1 ) GO TO 10 

I F ( I .EQ.IMAX-l ) GO TO 20 

****HERE THE NORMAL CONDITIONS ALONG THE N/S PERIMETER 
**** BETWEEN I = I M I N+2 AND I = IMAX-2 ARE GENERALLY STATED 

J= JM I N 
DELVX = 0 .0 

DELUY = (4*(U( I , J) +U( 1+1 , J) ) J -U { 1-1, J) -U ( I , J-l )-U( 1 + 1 
1 , J-l )-U( 1+2, J ) 

FI ( I , J)=F1 ( I, J )+{ ( DELVX- DELUY)/ (DEL TAX) )*A 
J= JMAX 
DELVX= 0 . 0 

DEL JY= ( -4* ( U ( I , J + 1 )+U( I + ljJ + 1) ) + {U( I-1»J + 1)+U( 1 + 2, J + l) 
1 + U( I , J + 2) + U< 1 + 1 , J+2) ) ) 

FI ( I , J ) = F1 { I , J ) + ( ( DELVX- DELUY) /(DELTAX) )*A 
GO TO 343 
10 J = J M I N 

DELVX = ( —V ( I— 1,J+1 ) —V { I— 1,J) J/2.0 

DELUY = ( U ( 1-1, J + l)-U( I - 1 , J ) )/2.0 + { 4* ( U ( I , J ) +U{ I + 1 , 

1J ) )-U( 1+2, J)-d( 1+1 , J-l ) -U( I , J-l ) ) 

Fl(I,J)=Fl(I,J)+(( DEL VX-DELUY) /(DELTAX) ) *A 
J= JMAX 

DELVX =(-V(I-l,J)-V(I-l,J+l) J/2.0 

DEL JY =(U( I — 1, J + l)— U( I — 1,J J ) /2.0 + { -4* l U ( I ,J+1)+U(I + 1 
1 , J + l ) )+U(I+2,J + l)+U(I + l,J+2)+U( I, J+2) ) 
Fl(I,J)=Fl(I,J)+{( CEL VX-DELUY) /(DELTAX) ) *A 
GO TO 343 
23 J = J M I N 

DELVX= ( V ( 1+2, J ) +V( 1+2, J+l) ) /2.0 

DELUY =(U(I+2»J + 1)-U(I+2,J) J/2.0 + ( 4* ( U ( I , J ) +U( 1 + 1 , J 
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I ) )-U( I + 1 , J-l)-U( I , J-l ) -U (1-1 , J) ) 
F1(I,J)=F1(I,J)+((DELVX-DELUY)/(DELTAX))*A 
J= JMAX 

DELVX = (V(I+2,J+1)+V(I+2,J) )/2.0 

DELUY = (U(I+2,J+1)-U(I+2,J) J/2.0 + ( -4* ( U (I , J + 1 ) +U< 

1 I + 1 , J + l ) )+U(I-l»J+l)+U ( I i J+2 ) +U ( I +1 , J+2 ) ) 
Fl(I,J)=Fl(I,J)+( ( DELVX-DELUY) /(DELTAX) )*A 
343 CONTINUE 

****THE FOLLOWING CALCULATIONS ARE FOR THE WESTERN ** 
****AND EASTERN PERIMETERS OF THE ISLAND , LESS CORNER 
**** NODAL POINTS . j*****^**#**#****#**#**:****;?##* 

****TEST FOR ZERO NODAL POINT ISLAND . **************** 

2 DO J J=J MAX— JM I N 

IF(JJ.LE.l) G3 TO 201 

MN=JMIN+1 

MX=J MAX-1 

DO 3 53 J = M N » M X 

IF(J.EQ.JMIN+1 ) GO TO 15 

IFIJ.EQ. JMAX-1 ) GO TO lfc 

**** HE p E THE FORMAL CONDITIONS ALONG THE W/E PERIMETER 
**** BETWEEN J = j M I N+2 AND J = J MAX- 2 ARE GENERALLY STATED 
I = I M I N 

DELVX = (4*( V( I , J ) +V( I , J + l) ) )-( V( I , J+2)+V( 1-1 , J+l) +V 
1(1-1 ,J)+V(I, J-l) ) 

DELUY=0. 0 

Fl(I,J)=Fl(I,J)+(( DELVX-DELUY) /(DELTAX) )*A 
I =1 MAX' 

DELVX = (-4*( V( I +1 , J+l )+V( 1 + 1 , J) ) )+V( 1+1 , J+2) +V( 1+2 , J+ 

II ) +V( 1+2, J )+V( 1+1, J-l) 

DE LUY=0 . 0 

F 1 ( I,J)=F1( I,J )+( (DELVX-DELUY) /(DELTAX) )*A 

GO TO 353 

15 I = I M I N 

DELVX = (V( 1 + 1, J-l ) -V ( I , J-l) ) /2.0 + (4*{ V( I , J)+V( I , J + l 
1 ) ) -V ( I , J+2)-V( I — 1 , J +1 ) —V ( I — 1 , J ) ) 

DELUY = (-U{ I , J-l) -U( I + 1 , J-l) ) /2.0 

FI ( I , J ) = F1( I , J ) + (( DELVX-DELUY )/(DELTAX) ) *A 

I = I MAX 

DELVX = (V( I + 1 , J-l )-V( I , J-l) ) /2.0 + ( -4* (V(I + 1,J)+V(I + 
1 1, J + l ) ) )+V ( I +1 , J+2 )+V( 1+2, J + l)+V( 1 + 2, J ) 

DELUY = C - UC 1+1 ,J-1)-U(I , J-l) J/2.0 

FI ( I, J )=F1( I , J )+( (DELVX-DELUY) /(DELTAX) )*A 

GO TO 353 

16 1= I M I N 

DELVX = (V( I+l,J+2)-V(I, J+2) J/2.0 + ( 4* ( V( I , J ) + V{ I , J+l 

I ) ) —V ( I ,J-1)-V( I — 1 , J ) —V ( I — 1,J+1 ) ) 

DELUY = (U( 1+1 , J+2) +U( I , J+2) ) /2.0 

Fl( I , J ) = FI ( I, J ) + ( ( DELVX-DELUY )/( DELTAX) ) *A 
I = I MAX 

DELVX = (V( I +1 , J+2 ) -V( I ,J+2) J/2.0 + ( -4* { V ( I +1 , J+ 1 ) 
1+V(I+1,J )) )+V( I + 1 , J-l ) +V ( I+2,J )+V( I + 2, J + l) 

DELUY = (U( I,J+2)+U(I+l,J+2) J/2.3 
F 1 ( I,J) = F1( I , J ) + ( ( DELVX-DELUY) /(DELTAX) )*A 
353 CONTINUE 

****COMMENCE THE ISLAND CORNER CALCULATIONS : ******** 
*v**CASE ONE ***** THE SOUTHWEST CORNER OF ISLAND. *** 

201 I = I M I N 
J = J M I N 

DELVX = (4*(V(I ,J)+V(I,J + 1) ) — (V (I — 1 ,J + 1 1 +V( I — 1,J)+V(I , J- 

I I) +V( I, J+2) ) J/2.0 

DE LUY= ( 4* ( U( I , J) +U ( 1 + 1 , J ) )-(U ( 1-1, J )+U ( I , J-l )+U( 1 + 1, J- 
1 1) +U ( I +2 , J ) ) ) / 2 . 0 

FI ( I , J ) = F1 ( I, J ) + ( ( DELVX-DELUY)/ (DELTAX) ) *A 

****CASE TWO ***** THE NORTHWEST CORNER OF ISLAND. *** 

I = I M I N 
J= JMAX 
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DEL VX= ( ( 4* (V(I,J)+V(I,J+l)))-(V(I»J+2)+V(I-l,J+l)+V(I- 
11,J)+V(I ,J-1) ) ) / 2 . 0 

DELUY= ( (-4*( U( I , J+l) + U( 1 + 1 , J + l ) ) )+U( 1-1 » J+1)+U(I , J+2) + 
1U( 1+1, J+2)+U( 1+2, J+l) )/2.0 
Fl(I,J)=Fl(I,J)+(( DELVX-OELUY )/ (DELT AX) )*A 

****CASE THREE *** THE NORTHEAST CORNER OF ISLAND. *** 

I = I MAX 
J= JMAX 

DE LVX= ( (-4*(V(I+1,J+1)+V(I+1,J ) ) )+V( 1 + 1, J + 2)+V( 1 + 2, J + l 
1)+V(I+2,J)+V(I+1,J-1) J/2.0 

DELUY= ( ( — A* ( U ( I,J+1)+U(I+1,J+1)))+U(I-1,J+1)+U(I , J+2J+ 
1U(I+1, J+2)+U( I +2, J+l) J/2.0 
FI ( I » J ) = F1 ( I » J ) + ( (DELVX-DELUY)/(DELTAX) )*A 

****CASE FOUR **** THE SOUTHEAST CORNER OF ISLAND. *** 

1= IMAX 
J= JMI N 

DELVX= ( (-4*( V( 1+1, J)+V( 1+1 , J+l) ) )+V( 1+1 , J + 2)+V(I+2, J + l 
1 )+V( 1+2 , J)+V( 1+1 , J-l) )/2.0 

DE LUY= ( ( 4* ( U( I ,J)+U(I+1,J)))-(U(I-1»J)+U(I»J-1)+U(I+1, 
1 J-l )+U( 1+2, J ) ) ) /2. 0 

FI (I »J)=F1 ( I , J ) + ( ( DELVX— DELUY)/(DELTAX) )*A 

************** AREA SIX CALCULATIONS: ***************** 
************** A REA $ j X IS THE AREA WHCIH DEPICTS **** 
************** ONLY THE INTERIOR OF THE ISLAND . ****** 
** TEST IS ESTABLISHED TO DETERMINE THE POSSIBILITY ** 
** OF A ZERO NODAL POINT ISLAND . ******************** 

I I=IMAX-IMIN 
J J=JMAX- JMI N 

IF( ( II.LE.D.OR. ( JJ.LE.l) ) GO TO 19 
NX= I MAX- 1 
N=I MI N+l 
DO 850 I =N , MX 
MN= J MIN+1 
MX=JMAX-1 
DO 850 J=MN , MX 
850 F1(I,J)=0.0 
19 RETURN 
END 



SUBROUTINE RELAX 



THIS SUBROUTINE PERFORMS THE RELAXATION PHASE OF THIS 
THESIS PROJECT . 

DVAR MEANS DIMENSIONED VARIABLES ************* 

COMMON/ ISLAND/ IM IN , IMAX, JMIN, JMAX, DPSIDI 
COMMON /DVAR/ PS I Ml , PS I , F 1 , U , V , D ELS QU , DELS QV , DPS I DT , 
1RESID, PSTEMP 

DIME NS ION PSIM1 (41 ,21 ) , PS I (41, 21 ) , FI (41 , 21 ) , U(42, 22) , 
1DELSQUI 40,20) , D E LS Q V( 40 , 20 ) , DPS IDT (41 ,21 ) , RES ID (41 , 21 ) 
2, PSTEMP (41, 21) ,V( 42,22) 

START THE RELAXATION 

SET CONSTANTS: EPSS AND ALPHA 

DE LT AX= 2. 0=* 10 . 0** ( 7) 

6=4.0*10. 0**8 
ALPHA=0 . 75 
WR IT E ( b , 105) ALPHA 
CF 1 =20 . 

EPSS= CFl*l.E-04 
WRITE(6, 106) EPSS 
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USE THE LAST KNOWN VALUE OF DPSIDT AS THE FIRST GUESS 



SET A SYSTEM OF ITERATIVE COUNTING. 

KC OUNT =0 
306 I COUNT = 0 

START THE COMPUTATION OF THE RESIDUAL. 

NOTE: DPS I DT= THE PARTIAL OF PSI WRT TIME, HAVING 

UNITS OF { CM. * *2 ) * (S EC .** {-2 ) ) . IT IS NOT THE 

PARTIAL OF PSI WRT TIME MULTIPLIED BY DELTAX SQUARED. 

NOTE: THE RESIDUAL IS THE LAPLACIAN OF DPSIDT 
MULTIPLIED BY THE SQUARE OF DELTAX, AND SUBTRACTING 
THE FORCING FUNCTION, FI . 

DO 444 JARE A=1 ,4 

I F ( { JMIN.EQ.2) . AND. (JAREA.EQ.l) ) GO TO 444 
IF( { JMAX.EQ.20) . AND. (JAREA.EQ.4) ) GO TO 444 
I F t { IMIN.EQ.2) . AND. (JAREA.EQ.2) ) GO TO 444 
IF{ { IMAX.EQ.40) . AN D . ( J AR EA . EQ . 3 ) ) GO TO 444 
IF( ( JAREA. EQ.I ) .OR. { JAREA. EQ.2 ) .OR. { JAREA. EQ.4) ) 

1 ILEFT= 2 

I F { (JAREA.EQ.l ) .OR . ( J AR E A . EQ . 3 } . OR . ( J AR EA . EQ . 4) ) 

HR I GHT =4 0 

I F { (JAREA.EQ.l ) ) JB0TM = 2 
L= JMI N-l 

IF(( JAREA.EQ.l J J JTOP=L 

I F( (JAREA.EQ.2 ) .CR . (JAREA.EQ.3 ) ) JBOTM=JMI N 
I F { (JAREA.EQ.2) . OR . { J A RE A . EQ . 3 ) ) JTOP=JMAX 
K= I M I N- 1 

IF ( (JAREA.EQ.2 ) ) I R IGHT = K 
KK=I MAX+ 1 

I F ( (JAREA.EQ.3) ) I L EFT = KK 
LL = JMAX + 1 

I F ( ( JAREA. EQ. A) ) JBOT M=LL 
I F { (JAREA.EQ.4) ) JT0P=20 
DO 400 J=JBOTM,JTOP 
DO 400 1=1 LEFT , I RIGHT 

RES I D{ I , J ) = { D^S IDT { 1+1 , J )+-DPSIDT{ I-1,J) + DPSIDT! I, J+l) 
1+DPSIDT (I , J-l )-4*DPS IDT( 1, J ) )-Fl( I , J ) 

CHECK FOR CONVERGENCE AND IF NECESSARY INITIALIZE 
VALUE OF CONSTANT C . 

THE REAL TEST FOR CONVERGENCE IS THE COMPARISON OF CF1 
AND RESID, THEREFORE' LE.EPSS ' MEANS 'MUCH LESS THAN 
CF1 • . 



I F ( ABS { RES ID( I , J ) ) .LE. EPSS ) GO TO 400 
C= UI.O + ALPHA) /4 .0 ) * (RESID( I,J ) ) 

ICOUNT = ICOUNT + 1 
DPSIDT (I, J)= DPSIDT! I , J )+C 
400 CONTINUE 

444 CONTINUE 

AREA FIVE CALCULATION . ************************** 

r5caftc * # # ape ^ * $ y. ^ ^ ^ 

DO 4A5 J=JMIN, JMAX 
DO 445 I =1 MIN, I MAX 

445 RESIDf I , J) = (DPSIDT( I + L, J) + DPSIDT{ 1-1, J) + DPSIDT( I , J + l) + 
1 DPSIDT! I , J-l )-4*DPS I DT ( I , J ) ) -FI ( I, J ) 

R E SB =0 . 0 
IMAXM1= I MAX— 1 
JMAX Ml = JMAX— 1 
RESI DI = 0. 0 



82 



ooooooooo oooooo 



DO 480 J=J M I N, JHAXM1 
DO 480 1=1 MI N, I MAX Ml 

RESB = 0.25*(RESID( I ,J)+RESID(I+1»J) +RE S I D ( I , J+l ) + 

1 RES I D ( 1+1, J + l) ) 

480 RESIDI =RESIDI + RE SB 

WRITE(6,107)KC0UNT,RESIDI,RESID(IMIN,JMIN),RESID(10,20 



CHECK FOR CONVERGENCE AND IF NECESSARY INITIALIZE 
VALUE OF CONSTANT Cl . 

I ICOUT = 0 

TERM = ( I MAX-I MIN + JMAX-JMIN) 

EPSI = E PS S 

IF(ABS(RESIDI) .LE.EPSI } GO TO 83 
ALPHAI =0.4 

Cl = (( 1.0+ALPHAI )/TERM)*(RESIDI) 

I COU NT = ICOUNT + 1 
I I COUT = 10 

DPSIDI = DPSIDI + Cl 
DO 490 J = J M I N, JMAX 
DO 490 I = 1 MI N» I MAX 
490 DPS I DT ( I , J ) = DPSIDI 
83 I F ( ICOUNT. EQ.O ) GO TO 500 
KCOUNT = KCOUNT + 1 
IFIKCOUNT.LT. 400) GO TO 306 
WRITE (6, 224) KCOUNT, ICOUNT, TKE, I ICOUT 
WRITE! 6,104) 

STOP ' 

PRINT OUT NEW DPSIDT ARRAY 



COMPUTATION THE TOTAL KINETIC ENERGY FOLLOWS: ******* 
THIS IS PLACED IN THE PROGRAM TO HELP INDICATE A ** 
STEADY-STATE CONDITION. **** 

500 DO 50 1=1,40 
DO 50 J=1 , 20 

AA= PSI (I+1,J+1) + PSI (I+1,J) 

BB= PSI ( I, J+l) + PSI (I,J) 

CC= PSI (I+1,J+1) +PSI (I, J+l) 

DD= PS I ( I + 1, J ) + PSI ( I , J ) 

EE = 2 . 0*DELT AX 

V ( I , J ) = ( (AA/EE) -( BB/EE) ) 

50 U(I,J) = ((DD/EE) -(CC/EE)) 

TKE=0. 0 
DO 100 1=1,40 
DO 100 J =1 , 20 

E = 0.5 * (U(I,J)**2 + V ( I , J ) **2 ) 

100 TKE=TK E + E 

WRITE (6 ,22 4) KCOUNT , ICOUNT, TKE, I ICOUT 
C WRITE! 6 , 102 ) 

DO 600 J = 2 , 2 0 
DO 600 1=2,40 

SCALE=10.0*=M-3) 

600 RESID! I , J )=DPS IDT! I ,J ) * SCALE 
C WR I TE ( 6 , 10 1 ) ( (RESID(I , J) ,1=1 ,33) ,J = 1,21) 

101 FORMAT! '0',///,21!3 3F4.2,//)) 

102 FORMAT (■ O' ,T63, ' RES ID F I EL D= DP S I DT FIELD * SCALE ',/// 
1 ) 

104 FORMAT (• ', T43, ' COMPUTATION TERMINATED DUE TO ', 

1' EXCESSIVE SWEEPS, — KCOUNT = 400 OR MORE ',//) 

105 FORMAT! 'O' , T63, ' AL P HA = ' ,F4.2,//) 

106 FORMAT ( ' 0 ' , T6 3 , ' EP S S= ' , F8 . 6 , // ) 

107 FORMAT! ' ' , I4.3F10.5) 

224 FORMAT! *0' , 'KCOUNT = ',14, 'ICOUNT = ',I4,'TKE= *,E10. 
14 , ' I ICOUT = ',14,//) 

RETURN 

END 

//GO. FT06F001 DD SYSOUT=A,SPACE=(CYL , ( 10, 2) ) 
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//GO .SYS IN DO * 

DATA CASE STUDY TWO-B USED. 

5 9 9 13 

J.V. SULLIVAN JR. STR EAMF UNC TI ON FIELD 
ISLAND-OCEAN MODEL ( I M IN= SIMAX= 9JMIN= 
J.V. SULLIVAN JR. STRE AMFUNCT I ON FIELD 
ISLAND-OCEAN MODEL (IMIN= 5IMAX= 9JMIN= 
J.V. SULLIVAN JR. ST RE AMFUNCT ION FIELD 
ISLAND-OCEAN MODEL! I MI N= 5IMAX= 9JMIN = 
J.V. SULLIVAN JR. ST RE AM FUNCT I ON FIELD 
ISLAND-OCEAN MODEL ( I M I N= 5IMAX= 9JMIN= 
J.V. SULLIVAN JR. STR EAMF UNCTION FIELD 
ISLAND-OCEAN MODEL! IMIN= 5IMAX= 9JMIN= 
J.V. SULLIVAN JR. STREAMFUNCT I ON FIELD 
ISLAND-OCEAN MODEL ( I M I N= 5IMAX= 9JMIN= 



PLOT 1 
9JMAX= 1 3 ) 
PLOT 2 
9 JMA X= 1 3 ) 
PLOT 3 
9 JMAX=1 3 ) 
PLOT 4 
9 JMAX= 13 ) 
PLOT 5 
9JMAX= 13) 
PLOT 6 
9 JMA X = 1 3 ) 
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